Lee- Yang zeros and phase transitions in nonequilibrium steady states 
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We consider how the Lee- Yang description of phase transitions in terms of partition function zeros 
applies to nonequilibrium systems. Here one does not have a partition function, instead we consider 
the zeros of a steady-state normalization factor in the complex plane of the transition rates. We 
obtain the exact distribution of zeros in the thermodynamic limit for a specific model, the boundary- 
driven asymmetric simple exclusion process. We show that the distributions of zeros at the first 
and second order nonequilibrium phase transitions of this model follow the patterns known in the 
Lee- Yang equilibrium theory. 
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With equilibrium statistical physics now firmly es- 
tablished as a successful general theory for predicting 
the macroscopic behavior of many-body systems, atten- 
tion has more recently turned to the the wider class of 
nonequilibrium systems. Examples of nonequilibrium 
systems are those relaxing to thermal equilibrium or 
those driven into a steady state, far from thermal equi- 
librium. A general framework for understanding these 
systems is proving elusive. However, over the last two 
decades the detailed study of specific models has revealed 
the wide range of phenomena that may emerge. 

It is by now well known that models of driven diffusive 
systems can exhibit phase transitions in their nonequilib- 
rium steady states Q. Typically such models comprise 
one or more species of particles being driven by an exter- 
nal field along with prescribed interactions between the 
boundaries and the outside world. Applications of these 
models are diverse and include the kinetics of biopoly- 
merization [|| , transport across a membrane || and traf- 
fic flow [|. From a more fundamental viewpoint, the 
interest lies in the richness of collective phenomena dis- 
played, including jamming H and spontaneous symmetry 
breaking even in one dimension. 

Half a century ago Lee and Yang provided a theory 
of equilibrium phase transitions based around the ze- 
ros of the partition function. In this work we focus on 
nonequilibrium steady states and how the phase transi- 
tions they admit may be placed in the context of the 
Lee- Yang theory. In the original exposition the fu- 
gacity of a gas was generalized to the complex plane and 
the zeros of the grand-canonical partition function in this 
plane were considered. In particular, it was shown that 
the free energy is analytic in any region of the complex- 
fugacity plane devoid of any zeros. Conversely, there are 
nonanalyticities at points where the density of partition 
function zeros accumulate in the thermodynamic limit. 
These accumulation points correspond to physically ob- 
servable phase transitions if they lie on the positive real 
axis. A similar phenomenon occurs in the complex plane 
of other intensive fugacity-like variables. For example, 



Lee and Yang also showed || that the zeros of the Ising 
model partition in the complex plane of activity (defined 
as exp(—h/T) where h is magnetic field and T tempera- 
ture) lie on the unit circle. Similarly, partition function 
zeros in the complex temperature plane also accumulate 
at phase transition points || [l0|. The applicability of 
the Lee- Yang theory to equilibrium transitions driven by 
intensive field-like variables can be understood from the 
fact that, mathematically, these variables play similar 
roles in the partition function. 

The difficulty in applying the Lee- Yang theory to, e.g., 
the steady state of a driven diffusive system is that one 
does not have to hand a partition function defined in 
terms of thermodynamic state variables. However, one 
can always define a quantity Z as a sum of the steady 
state weights (unnormalizcd probabilities) /(C): 
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Thus the probability of a configuration is given by 
P(C) = f(C)/Z. In this work, we treat the normalization 
Z as a nonequilibrium analog of the partition function. 

In order to obtain Z one still has to calculate the steady 
state weights. Up to a multiplicative factor they are im- 
plied by the transition rates W(C — > C) that define the 
model through the requirement that the total inflow of 
probability into a given configuration C must be balanced 
by the outflow, i.e., 
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C) - f(C)W(C -► C')} = . (2) 



The solution of this set of equations for the steady-state 
weights f{C) can, in principle, be obtained using stan- 
dard methods, such as Gaussian elimination, Cramer's 
rule or graphically If one applies such a method, 

the weights can always be written as polynomials of the 
elementary transition rates. Thus, for a finite system, Z 
is also a polynomial of the transition rates. 

In this work we generalize the transition rates to the 
complex plane and in doing so consider the zeros of the 
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FIG. 1: Dynamics of the ASEP. The arrow labels indicate the 
rates at which the corresponding transitions occur; site labels 
are also indicated. 



partition function. In analogy with the Lee- Yang picture 
for equilibrium systems, for a phase transition to arise we 
expect the zeros to pinch the real axis at a point which 
corresponds to a (real) critical transition rate. 

As an aside, we note that the general methods of solv- 
ing (J2J) mentioned above are rarely tractable when the 
number of configurations becomes large. One exception 
is when the system satisfies detailed balance. In this case, 
each term within the square brackets in (||) is identically 
zero and it is easy to construct the set of steady-state 
weights as Boltzmann weights with respect to an energy 
function. The resulting Z is then an equilibrium parti- 
tion function. A simple example of such a system is the 
Ising model evolving under Glauber dynamics wherein 
studying the zeros of the partition function as a function 
of the transition rate would be equivalent to studying the 
complex-temperature zeros. 

In the present work, however, we are interested in mod- 
els which do not satisfy detailed balance and hence for 
which the normalization Z is not a priori known. This 
makes it difficult to make general statements about the 
significance of the normalization Z. We have argued that 
a nonequilibrium phase transition will be signaled by an 
accumulation of the zeros of Z in the complex plane of 
transition rates towards the real axis. We demonstrate 
this in the case of a particular nonequilibrium model for 
which an expression for Z is exactly known and we can 
calcuate the distribution of zeros analytically. We find 
that these zeros lie on circles in the complex plane of an 
appropriately chosen function of a transition rate that 
plays the role of a fugacity. The characteristic properties 
of the distributions of partition function zeros at first 
and second order equilibrium transitions |?], [h], |T^| are 
also exhibited in this nonequilibrium case. In short, we 
show how the the Lee- Yang theory of phase transitions 
generalizes to a nonequilibrium steady state. 

Before discussing the model and our results in greater 
detail, we wish to remark how our study differs from 
a related work due to Arndt [jl3|. The aim of that in- 
vestigation was to locate a phase transition induced by 
varying the relative numbers of different particle species 
in a particle-conserving driven diffusive system. To this 
end, a fugacity-like quantity was introduced in an ad hoc 
way. Physically, this corresponds to placing the driven 
diffusive system in contact with a particle reservoir thus 
introducing equilibrium fluctuations in the particle den- 
sities. We stress that in the present work, it is elementary 
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FIG. 2: Phase diagram of the ASEP. 



transition rates, which need not have any connection to 
intensive field (fugacity-like) variables, that are general- 
ized to the complex plane. 

We now introduce the model we consider, namely 
the one-dimensional asymmetric simple exclusion process 
(ASEP) with open boundaries. Since its introduction 
Q this model has received a great deal of attention — 
for further details and references. 



see, e.g., 11, |15| for further details and references. The 
model comprises a lattice of N sites, each of which may 
be occupied or empty. In an infinitesimal time interval dt 
one of the following transitions may occur: a particle may 
hop one site to the right with probability dt (i.e. at unit 
rate) subject to the receiving site being vacant; a parti- 
cle may be inserted onto the leftmost site (if vacant) with 
probability adt or a particle occupying the rightmost site 
(if present) may be removed with probability /3dt. Fig. |l| 
illustrates the microscopic dynamics. 

As outlined above, there is a current (defined as the 
flux of particles between neighbouring sites per unit time) 
in the steady state which in the thermodynamic limit 
N — > 00 exhibits nonanalyticities as the boundary rates 
a and (3 are varied. These are associated with phase 
transitions, and the phase diagram is given in Fig. ||. 

Phase (i) is a high-density phase with a density pro- 
file that decays exponentially towards the right bound- 
ary. In this phase the current J = [3{1 — 0). Since the 
ASEP is invariant under particle-hole exchange and the 
swap a <-> j3, it follows that phase (ii) is a low-density 
phase that has an exponential decay in the density pro- 
file from the left boundary and a current J = a(l — a). 
At the first order transition line a — (3 < | the current 
exhibits a discontinuity in its first derivative. Here the 
system exhibits a shock front separating regions of high 
and low density; this is an example of phase coexistence 
at a nonequilibrium phase transition. 

In phase (iii) , the current assumes a constant value J = 
2 which is the largest possible current admitted for any 
combination of a and /3. Hence phase (iii) is called the 
maximal current phase. This phase has the interesting 
feature that the density profile decays as a power-law 
from both boundaries and hence at the transition from 
either phase (i) or (ii) to phase (iii) a correlation length 
diverges. Additionally, J has a discontinuity in its second 



FIG. 3: Zeros of the normalization (^|) in the complex-a plane FIG. 4: Zeros of the normalization (^) in the complex-^ 
for j3 = |, 1. In both cases the lattice size N is 300. a(l — a) plane for (3 = |, 1 and N = 80, 160, 320. 



derivative and so this is a second-order transition. 

These results are known mainly from the exact solution 
Jl6| , [I?]]. The two key results we call upon here are that 
the normalization for an iV-site system Zn is given by 
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and the current by J = Z^-i/Z^- In order later to con- 
struct the distribution of the zeros of Zm for the ASEP 
we also need the large- TV forms of Zm- From the exact 
solution it can be shown that for all a, /3, Zjv ~ AJ~ N N 1 
where J is the current and A, 7 depend only on the phase 
being considered. We see In J is the extensive part of 
— In Zn and thus In J plays the role of the free energy 
and J the fugacity for this nonequilibrium system. 

We now consider the zeros of Zn given by (||) in the 
complex-a plane at fixed (3 18 1. (The symmetry in a 
and /3 of Zjq implies that one would obtain the same 
pattern of zeros in the complex-/? plane at fixed a) . Using 
Mathematica, we obtained numerical estimates of the 
zeros of Zn for system size (and number of zeros) N = 
300 and (3 = 1, \. The results are shown in Fig. H. 

For the case (3 = -|, the zeros appear to approach the 
positive real a axis at a value a c — | , coincident with the 
first-order phase transition point in this regime. The plot 
for the case (3 = 1 is suggestive of a slow accumulation 
of the zeros to the positive real a axis at a c = ^, i.e. the 
second-order phase transition point. It appears also in 
this case that the zeros approach the real axis at an angle 
j as opposed to an angle ^ at the first-order transition. 
We will shortly demonstrate that the above assertions — 
consistent with the patterns of zeros at first- and second- 
order phase transitions known from the equilibrium Lee- 
Yang theory [0, [l0| [l2| — are correct. 

It is useful to make a change of variable from a to £ = 
a(\ — a) (note that in the low density phase £ = J) and 
we plot the numerically-obtained zeros in the complex-^ 
plane in Fig. [|. For both (3 ~ it appears that, as 
N — > 00, the zeros become uniformly distributed around 
a circle of radius £ c = a c (l — a c ) (these circles are shown 



as solid lines in the figure). We shall now show how to 
obtain this result analytically. 



We first note (for a proof see |19 ) that for any poly- 
nomial pat(£) with N zeros, its density of zeros p(x,y) 
in the complex £ = x + iy plane is given by p(x, y) = 
2ttV 2 In \pn(0\ where V 2 = + -j^. In the thermody- 
namic limit (N — > 00) we expect a continuous distribu- 
tion of zeros and wish to normalize p so that its integral 
over all space is unity. We thus define the key quantity 
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so that p — 27rV 2 0. In equilibrium statistical physics 
</> would be interpreted as the extensive part of the free 
energy; as noted above cf) = — In | J| for the ASEP. 

It is useful to regard cf> as an 'electrostatic' potential 
in a two-dimensional space. We introduce the field E = 
V</> and consider the discontinuities that arise as a line 
of charge (zeros) is crossed. The appropriate forms of 
Gauss's and Ampere's laws imply that the component 
of E parallel to a line of zeros separating two phases 
must be continuous across it whereas the perpendicular 
component differs by a value 2-rtg{x,y) where g{x,y) is 
the line density of charge (zeros) along the boundary. 

The current itself has only two analytic forms. In the 
interior region of Fig. [I] (i.e., the region that includes 
£ = 0) J = a(l — a) = £ and hence <f> = — In |£| . Over the 
exterior region, the current is constant: J = a c (l — a c ) — 
£ c and so tfi = — ln|£ c |. 

The field that results from the potential <f> is zero in the 
exterior region and directed radially towards the origin 
in the interior region; specifically E = — f/r 2 where r = 
(Re£,Im£). Now, since the component of the electric 
field parallel to the boundary between the two phases 
varies continuously and is zero in the exterior region, the 
phase boundary must lie perpendicular to the field lines 
in the interior region. In other words, the phase boundary 
is an equipotential ^ = — ln.|£| = const =4> |£| = £ c since 
the boundary must pass through the transition point at 
£ = £ c . The line-density of zeros g(0) along the circle £ = 
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£, c e l9 can be calculated by noting that the perpendicular 
component of the field at the boundary has a magnitude 
l/|£ c |- The relationship given above then implies g(9) — 
(27r£ c ) _1 , i.e., a constant as claimed. 

We now use this result to investigate the distribution 
of zeros near the transition point in the complex-a plane. 
In the £ plane, the circle of zeros passes along the line 
Rc £ = £ c which maps onto a curve in the complex a — 
u + iv plane that has (u ~ \) 2 — v 2 = i — £ c . The solution 
of this equation that describes a line passing through the 
transition point a = a c is u = | — (v 2 + \ — ^c) 1 ^ 2 - 

When the phase transition is first order, (3 < \ and 
£ c = (3(1 — (3) < j. Then, the curve of zeros passes 
smoothly through the transition point a = f3. One can 
show that the line density of zeros at this point in the 
complex-a plane is (1 — 2(3) /[2tt(3(1 — (3)], i.e., nonzero. 
The fact that the density of zeros is nonzero on the real 
axis at a first-order phase transition is well-known in the 
equilibrium Lee- Yang theory [Q. 

At the second order phase transition ((3 > i and £ c = 
j) the line of zeros is given by u = h — \v\. This means 
the zeros approach the transition point a = 1/2 at an 
angle of ? to the real axis, meeting at a right angle. If 
one defines I as the displacement along the straight line 
of zeros from the transition point, one finds that it varies 
with y — Im£ as £ = y/y. This implies that the density 
of zeros behaves as g(£) — ~4| = ^r. In the equilibrium 
theory, it is well known that the density of zeros vanishes 
as a power-law towards the transition point on the real 
line (see, e.g., Q), a result we have recovered in this 
nonequilibrium case. 

To summarize, we have shown that the zeros for the 
normalization of the ASEP accumulate at the phase tran- 
sition points in the complex plane of transition rates. As 
in the case of the Lee- Yang theory, a first-order transi- 
tion is manifested by a nonzero density of zeros at the 
accumulation point, whereas the density of zeros decays 
as a power-law towards a continuous transition point. 

A point that we have so far ignored is that the normal- 
ization (Q) is only defined up to a multiplicative factor 
that depends on the method used to solve Eqn. (||). This 
additional factor is analogous to that which would ap- 
pear in an equilibrium partition function after a uniform 
shift of the energy scale. Although this factor could itself 
be a polynomial of the transition rates and so introduce 
additional zeros, one would not expect these spurious ze- 
ros to be physically relevant to the phase behavior. For 
example one way to define the normalization is to use 
Cramer's rule to solve (§) (see JO)). Then for the ASEP 
the normalization is given by (a/3) N times the expression 
and additional zeros are introduced only at the origin. 
Finally we remark on the applicability of the Lee- Yang 
theory to other nonequilibrium systems. Firstly, a gen- 
eralization of the ASEP that allows particles also to hop 
to the left at a rate q has been solved J2(J|. The phase 
behavior when q < 1 is very similar to that of the ASEP, 



and the rate q enters into the expressions for the current 
in a simple way. The result of this would be for the pat- 
tern of zeros in the a-plane to shrink towards the origin 
as q — ► 1. At this point, there is a transition to a regime 
where the current ceases to flow in the thermodynamic 
limit. It would be of interest to see whether the zeros 
of the normalization in the complex-g plane accumulate 
at q = 1 since this would provide some evidence for the 
generality of the Lee- Yang picture when applied to com- 
plex transition rates. Further evidence is provided by 
a study of the percolation probability on finite directed 
percolation lattices |2lJ] in which a similar accumulation 
phenomenon to that reported here was observed. Also, 
one might learn something of the approach to the ther- 
modynamic limit by considering how the distribution of 
zeros varies as the system size is increased. 
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